# Current, as of Thursday, May 28, 2020.
# rm(list = ls())
# setwd("~/Dropbox/Base Dropbox/List Experiment Paper (PS)/Conditional Accept/Replication Files (Carmines-Schmidt PS)/2012 CCAP/Partisanship & Ideology Models")
library(foreign)
library(car)
library(readstata13)
library(zeligverse)
library(ggplot2)
#-----
# Input the data.
our.data = read.dta13("CarminesSchmidtPS-replication-2012.dta", generate.factors = TRUE)
names(our.data)
#-----
# Subset the data into only respondents that were assigned to either the BASELINE or 'MORMON' treatment conditions. 
mormons.experiment = subset(our.data, treatment=="Mormons" | treatment=="Control")
#-----
# Model: treatment effects as function of PARTISANSHIP. 
our.model = lm(item_count ~ mormon_binary + republican + democrat + 
                 mormon_binary*republican + mormon_binary*democrat, 
               data = mormons.experiment)
summary(our.model)

# Simulate: treatment effects for REPUBLICANS.
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$republican*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$democrat*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for REPUBLICANS.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + republican + democrat + 
                interaction1 + interaction2,  
              model = "normal", data = mormons.experiment)
z.out
# REPUBLICANS in the BASELINE condition. 
x.low = setx(z.out, 
             mormon_binary = 0,
             republican = 1,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0)
# REPUBLICANS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              republican = 1,
              democrat = 0,
              interaction1 = 1,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for REPUBLICANS. 
LB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for DEMOCRATS. 
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$republican*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$democrat*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for DEMOCRATS. 
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + republican + democrat +
                interaction1 + interaction2,  
              model = "normal", data = mormons.experiment)
z.out
# DEMOCRATS in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             republican = 0,
             democrat = 1,
             interaction1 = 0,
             interaction2 = 0)
# DEMOCRATS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              republican = 0,
              democrat = 1,
              interaction1 = 0,
              interaction2 = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for DEMOCRATS. 
LB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for INDEPENDENTS.
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$republican*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$democrat*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for INDEPENDENTS.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + republican + democrat +
               interaction1 + interaction2,  
              model = "normal", data = mormons.experiment)
z.out
# INDEPENDENTS in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             republican = 0,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0)
# INDEPENDENTS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              republican = 0,
              democrat = 0,
              interaction1 = 0,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for INDEPENDENTS. 
LB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Model: treatment effects as function of IDEOLOGY.
our.model = lm(item_count ~ mormon_binary + conservative + liberal +
                 mormon_binary*conservative + mormon_binary*liberal, 
               data = mormons.experiment)
summary(our.model)

# Simulate: treatment effects for CONSERVATIVES.  
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$conservative*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$liberal*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for CONSERVATIVES in the treatment condition. 
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + conservative + liberal + 
                interaction1 + interaction2,  
              model = "normal", data = mormons.experiment)
z.out
# CONSERVATIVES in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             conservative = 1,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0)
# CONSERVATIVES in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              conservative = 1,
              liberal = 0,
              interaction1 = 1,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for CONSERVATIVES. 
LB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for LIBERALS.
# Set up the interaction terms. 
interaction1 = mormons.experiment$conservative*mormons.experiment$mormon_binary
interaction2 = mormons.experiment$liberal*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for LIBERALS.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + conservative + liberal +
                interaction1 + interaction2,  
              model = "normal", data = mormons.experiment)
z.out
# LIBERALS in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             conservative = 0,
             liberal = 1,
             interaction1 = 0,
             interaction2 = 0)
# LIBERALS in the MORMON condition. 
x.high = setx(z.out, 
              mormon_binary = 1,
              conservative = 0,
              liberal = 1,
              interaction1 = 0,
              interaction2 = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for LIBERALS.   
LB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for MODERATES.
# Set up the interaction terms. 
interaction1 = mormons.experiment$conservative*mormons.experiment$mormon_binary
interaction2 = mormons.experiment$liberal*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for MODERATES.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + conservative + liberal +
                interaction1 + interaction2,  
              model = "normal", data = mormons.experiment)
z.out
# MODERATES in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             conservative = 0,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0)
# MODERATES in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              conservative = 0,
              liberal = 0,
              interaction1 = 0,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for MODERATES. 
LB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Point estimates.
# round(cbind(PE.lib, PE.mod, PE.con, PE.dem, PE.ind, PE.gop), digits = 2)
# Lower bounds.
# round(cbind(LB.lib, LB.mod, LB.con, LB.dem, LB.ind, LB.gop), digits = 2)
# Upper bounds.
# round(cbind(UB.lib, UB.mod, UB.con, UB.dem, UB.ind, UB.gop), digits = 2)
#-----
# Model: treatment effects as function of PARTISANSHIP and additional covariates. 
our.model = lm(item_count ~ mormon_binary + republican + democrat +
                 republican*mormon_binary + democrat*mormon_binary +
                 ideology + female + white + black + education + church + rr +
                 obama_favor + interest + income + romney_favor,  
               data = mormons.experiment)
summary(our.model)

# Simulate: treatment effects for REPUBLICANS, controlling for confounding variables. 
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$republican*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$democrat*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for REPUBLICANS.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + republican + democrat +
                interaction1 + interaction2 +
                ideology + female + white + black + education + church + rr +
                obama_favor + interest + income + romney_favor,  
              model = "normal", data = mormons.experiment)
z.out
# REPUBLICANS in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             republican = 1,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0,
             ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# REPUBLICANS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              republican = 1,
              democrat = 0,
              interaction1 = 1,
              interaction2 = 0,
              ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for REPUBLICANS. 
LB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate, treatment effects for DEMOCRATS, controlling for confounding variables. 
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$republican*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$democrat*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for DEMOCRATS. 
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + republican + democrat +
                interaction1 + interaction2 +
                ideology + female + white + black + education + church + rr +
                obama_favor + interest + income + romney_favor,  
              model = "normal", data = mormons.experiment)
z.out
# DEMOCRATS in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             republican = 0,
             democrat = 1,
             interaction1 = 0,
             interaction2 = 0,
             ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# DEMOCRATS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              republican = 0,
              democrat = 1,
              interaction1 = 0,
              interaction2 = 1,
              ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for DEMOCRATS. 
LB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate, treatment effects for INDEPENDENTS, controlling for confounding variables. 
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$republican*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$democrat*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for INDEPENDENTS.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + republican + democrat +
                interaction1 + interaction2 +
                ideology + female + white + black + education + church + rr +
                obama_favor + interest + income + romney_favor,  
              model = "normal", data = mormons.experiment)
z.out
# INDEPENDENTS in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             republican = 0,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0,
             ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# INDEPENDENTS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              republican = 0,
              democrat = 0,
              interaction1 = 0,
              interaction2 = 0,
              ideology = mean(mormons.experiment$ideology, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for INDEPENDENTS. 
LB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Model: treatment effects as function of IDEOLOGY and additional covariates. 
our.model = lm(item_count ~ mormon_binary + conservative + liberal +
                 mormon_binary*conservative + mormon_binary*liberal +
                 pid_7 + female + white + black + education + church + rr +
                 obama_favor + interest + income + romney_favor, 
               data = mormons.experiment)
summary(our.model)

# Simulate: treatment effects for CONSERVATIVES, controlling for confounding variables. 
# Set up the interaction terms. 
mormons.experiment$interaction1 = mormons.experiment$conservative*mormons.experiment$mormon_binary
mormons.experiment$interaction2 = mormons.experiment$liberal*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for CONSERVATIVES. 
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + conservative + liberal + 
                interaction1 + interaction2 + pid_7 + female + white + black +
                education + church + rr + obama_favor + interest + income + 
                romney_favor,
              model = "normal", data = mormons.experiment)
z.out
# CONSERVATIVES in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             conservative = 1,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0,
             pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# CONSERVATIVES in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              conservative = 1,
              liberal = 0,
              interaction1 = 1,
              interaction2 = 0,
              pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for CONSERVATIVES.
LB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for LIBERALS, controlling for confounding variables.
# Set up the interaction terms. 
interaction1 = mormons.experiment$conservative*mormons.experiment$mormon_binary
interaction2 = mormons.experiment$liberal*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for LIBERALS.
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + conservative + liberal +  
                interaction1 + interaction2 + pid_7 + female + white + black +
                education + church + rr + obama_favor + interest + income + 
                romney_favor,
              model = "normal", data = mormons.experiment)
z.out
# LIBERALS in the BASELINE condition. 
x.low = setx(z.out, 
             mormon_binary = 0,
             conservative = 0,
             liberal = 1,
             interaction1 = 0,
             interaction2 = 0,
             pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# LIBERALS in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              conservative = 0,
              liberal = 1,
              interaction1 = 0,
              interaction2 = 1,
              pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for LIBERALS. 
LB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for MODERATES, controlling for confounding variables. 
# Set up the interaction terms. 
interaction1 = mormons.experiment$conservative*mormons.experiment$mormon_binary
interaction2 = mormons.experiment$liberal*mormons.experiment$mormon_binary
# Use Zelig to model the difference-in-means for MODERATES. 
set.seed(47408)
z.out = zelig(item_count ~ mormon_binary + conservative + liberal + 
                interaction1 + interaction2 + 
                pid_7 + female + white + black + education + church + rr +
                obama_favor + interest + income + romney_favor,
              model = "normal", data = mormons.experiment)
z.out
# MODERATES in the BASELINE condition.
x.low = setx(z.out, 
             mormon_binary = 0,
             conservative = 0,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0,
             pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
             female = median(mormons.experiment$female, na.rm=TRUE),
             white = median(mormons.experiment$white, na.rm=TRUE),
             black = median(mormons.experiment$black, na.rm=TRUE),
             education = median(mormons.experiment$education, na.rm=TRUE),
             church = median(mormons.experiment$church, na.rm=TRUE),
             rr = mean(mormons.experiment$rr, na.rm=TRUE),
             obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
             interest = median(mormons.experiment$interest, na.rm=TRUE),
             income = median(mormons.experiment$income, na.rm=TRUE),
             romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# MODERATES in the MORMON condition.
x.high = setx(z.out, 
              mormon_binary = 1,
              conservative = 0,
              liberal = 0,
              interaction1 = 0,
              interaction2 = 0,
              pid_7 = mean(mormons.experiment$pid_7, na.rm=TRUE),
              female = median(mormons.experiment$female, na.rm=TRUE),
              white = median(mormons.experiment$white, na.rm=TRUE),
              black = median(mormons.experiment$black, na.rm=TRUE),
              education = median(mormons.experiment$education, na.rm=TRUE),
              church = median(mormons.experiment$church, na.rm=TRUE),
              rr = mean(mormons.experiment$rr, na.rm=TRUE),
              obama_favor = median(mormons.experiment$obama_favor, na.rm=TRUE),
              interest = median(mormons.experiment$interest, na.rm=TRUE),
              income = median(mormons.experiment$income, na.rm=TRUE),
              romney_favor = median(mormons.experiment$romney_favor, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for MODERATES.
LB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Point estimates for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(PE.lib, PE.mod, PE.con, PE.dem, PE.ind, PE.gop), digits = 2)
# Lower bounds for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(LB.lib, LB.mod, LB.con, LB.dem, LB.ind, LB.gop), digits = 2)
# Upper bounds for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(UB.lib, UB.mod, UB.con, UB.dem, UB.ind, UB.gop), digits = 2)



